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Abstract 

We study a three-parameter family of "PT-symmetric Hamiltonians, related via the ODE/IM 
correspondence to the Pcrk-Schultz models. We show that real eigenvalues merge and become 
complex at quadratic and cubic exceptional points, and explore the corresponding Jordon 
block structures by exploiting the quasi-exact solvability of a subset of the models. The 
mapping of the phase diagram is completed using a combination of numerical, analytical 
and perturbative approaches. Among other things this reveals some novel properties of the 
Bender-Dunne polynomials, and gives a new insight into a phase transition to infinitely-many 
complex eigenvalues that was first observed by Bender and Boettcher. A new exactly-solvable 
limit, the inhomogeneous complex square well, is also identified. 
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1 Introduction 



In this paper we return to the spectra of a family of "PT-symmetric eigenvalue problems first 
studied in detail in [1,2]. Consider the following differential operator: 
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where M, a and A are real numbers, with M > 0, and the powers of ix are rendered single- 
valued by placing a cut along the positive imaginary x axis. Then an eigenvalue problem, with 
a discrete spectrum, can be defined as 



Htp(x) = E^{x); ip(x)eL 2 (C), 



;i.2) 



where C is an infinite contour in the complex plane, which must pass below the origin whenever 
A 2 ^ j or M ^ Z. For M < 2 the ends of this contour can asymptote to the negative and 
positive real axes, while for M > 2 they must be deformed down into the complex plane so as 
to continue the M < 2 spectral problem smoothly [3]. This is illustrated in figured] 
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Figure [T] A possible quantisation contour C for M just 
larger than 2, together with some of the Stokes sectors. 

In [1, 2], the eigenvalue problem was initially specified in terms of I = A — i, but with 
boundary conditions imposed at infinity, the choice of A is more natural. 

An alternative specification of the boundary conditions, which holds for all values of M, 
starts from the Stokes sectors for (jl.ip . which we denote by 
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For all M the requirement (|1.2[) is equivalent to the demand that ip(x) — > as x — > oo in the 
sectors 5_i and Si. This allows for a convenient rephrasing of the eigenvalue condition in terms 
of the vanishing of a certain Wronskian. Following Hsieh and Sibuya [4-6], let y (x,E,a,X) be 
the solution to (JTTTJ) that is (uniquely) defined by the following asymptotic for x — > oo on the 
negative imaginary axis: 



V (x, E, a, A) ~ ± {ix)- M '*-^ exp (-^^ 



x — > — l oo 



and then set 



u = e 



i7r/(M+l) 



(1.4) 
(1.5) 
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and define a sequence of further solutions y^ to (jl.ip by 

y k (x, E, a, A) = ^/2-(-i)*W2 y Q ( u -k X] ^^E, (-l) fc a, A) . (1.6) 

It is easily checked that y^ decays, or is subdominant, in Su, and that the 'nearest-neighbour' 
Wronskians W[yk,yk+i] = VkVk+i ~ y'kVk+i are an equal to 10 The eigenvalue condition is then 
that y_i and y\ should be proportional to each other, in other words that E should be a zero 
of the 'next-nearest-neighbour' Wronskian 

T(E,a,\) = W[y- 1 ,y 1 ). (1.7) 

From this characterisation, and the analyticity of T as a function of its arguments, a number of 
important properties, such as the discreteness of the spectrum, immediately follow. In addition 
to being a spectral determinant, via the ODE/IM correspondence of [7] (see [8] for a review) T 
encodes the properties of the ground state of an integrable quantum field theory, in this case the 
Perk-Schultz model [9,10]. This correspondence is based in part on the fact that T is a Stokes 
multiplier for (jl.ip . in that the following equation holds [1,11]: 

T(E, a, X)y (x, E, a, A) = y_i(x, E, a, A) + y 1 (x, E, a, A) . (1.8) 

A feature the eigenvalue problem (II. 2h shares with many other 'PT-symmetric problems is 
the reality of its spectrum for many values of the free parameters. In particular, for real M > 1, 
a and A, the spectrum of (|l.lj) can be proved to be 

• real if a < M + 1 + 2 |A| ; (1.9) 

• positive if a < M + 1 - 2 |A| . (1.10) 

These results were established in [1] using techniques inspired by the ODE/IM correspondence. 
One of the main aims of this paper is to refine this picture and to explore in more detail how 
and where spectral reality is lost as the region (jl.9p is left. 

Along the lines a = M + 1 ± 2A which form the frontiers of the region (|1.9j) of guaranteed 
reality, the model has an exactly-zero energy level, as in supersymmetric quantum mechanics. 
The 'protection' of this level can be seen as the mechanism by which the first levels become 
complex [2]. However, numerical investigations at M = 3, reported in [2], showed that the region 
within which the spectrum of (jl.ip is complex has considerably more structure than (jl.9p might 
suggest. The curved, cusped line of figure [2] indicates where the first pair of complex eigenvalues 
is formed as the region of complete spectral reality is left; it touches the lines a = M + 1 ± 2A 
at isolated points, where the protected zero-energy level coincides with another level. 

*Note, a propagating typo in [1] and [8] resulted in the factor of (— l) fc multiplying ka/2 in the exponent of uj 
being omitted from the definition of yt given in those papers. None of the other formulae in [1,8] are affected. 
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Figure |21 The domain of unreality in the (2 A, a) plane for M = 
3, with portions of lines with a protected zero-energy level also 
shown. The horizontal axis is 2A = 21 + 1. 

The additional dotted lines on the figure, also at angles of ±45°, show points within the 
region a > M + 1 + 2|A| where the model has an exactly- zero energy level; exceptionally for 
M = 3, the model is also quasi-exactly solvable along these lines. It is notable that, to within 
numerical accuracy, the cusps on the boundary of the region of unreality appear to lie exactly 
on these lines. 

It is natural to ask where further pairs of complex eigenvalues are formed. For M = 3 the 
answer is shown in figure El adapted from [8]; the same pattern was found independently by 
Sorrell [12] via a complex WKB treatment of the problem. The pattern of cusps is repeated, 
with the cusps again appearing to lie on the lines of protected zero-energy levels. 
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Figure [3] The (2A, a) plane for M = 3, showing lines across 
which further pairs of complex eigenvalues are formed. 
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The analysis of [2] left a number of questions open. Whilst the merging and subsequent 
complexification of levels was suggestive of exceptional points and a Jordan block structure for 
the Hamiltonian, this was not demonstrated explicitly. The apparent siting of the cusps on lines 
with simultaneous quasi-exact solvability and protected zero-energy levels was not proved; in 
particular, it was not clear whether this feature should be associated with the zero-energy level 
(in which case it should persist for M / 3) or with the quasi-exact solvability (in which case it 
might be lost for Finally, and connected with this last question, the general pattern 

away from M = 3 was not explored. 

In this paper we revisit these issues. For M = 3 we investigate the positions of the cusps, 
proving that they do indeed lie on QES lines, and look at the exceptional points in the spectrum 
and the Jordan form at such points. We then explore the situation for M^3 numerically, and 
verify the picture that emerges with detailed perturbative studies near M = 1 and M = oo. The 
perturbative treatment near M = 1 also gives a new insight into the transition to infinitely-many 
complex eigenvalues for M < 1, first observed by Bender and Boettcher for the a = 0, A 2 = \ 
case of (jl.ip . 

2 Exact locations of special exceptional points 

2.1 Generalities and previous results 

Exact formulae for the full curves of exceptional points are unlikely to exist, even for M = 3. 
However, certain exceptional points can be located exactly, and this information turns out to be 
very useful in mapping the full phase diagram. As in [2], we begin by introducing an alternative 
set of coordinates on the (2A, a) plane, defined by 



For M = 3 these coordinates are illustrated in figureUl The lines a+ G N and a_ G N correspond 
to the dotted lines on figures [2] and [3l along which the model (jl.ip has an exactly-zero energy 
level. For M = 3 the existence of this level can be understood in terms of quasi-exact solvability 
and a hidden AA-fold supersymmetry [1, 13]. There are also exactly-zero energy levels along the 
lines a± = 0, related for all values of M to standard quantum-mechanical supersymmetry [2]. 

Exceptional points occur in the spectrum of an eigenvalue problem whenever the coalescence 
of two or more eigenvalues is accompanied by a coalescence of the corresponding eigenvectors; 
at such points there is a branching of the spectral surface [14-17]. In 'PT-symmetric systems, 
eigenvalues are all either real, or in complex-conjugate pairs. Complex eigenvalues can therefore 
be formed only via the intermediate coincidence of two (or more) previously-real eigenvalues. 
For one-dimensional problems of the sort under discussion here genuine degeneracies of levels 
are impossible - since, for example, the Wronskian of any two solutions which both decay 
exponentially in the same asymptotic direction must vanish - and so levels in our problem can 
only coincide at exceptional points. Hence the cusped lines on figures [2] and [3] are lines of 
exceptional points. In fact, we shall see that points on the (codimension one) smooth segments 
of the cusped lines are quadratically exceptional, with two levels coalescing, while the cusps 
themselves, of codimension two, are cubic exceptional points. In the following, we will often 
refer to a connected union of quadratically-exceptional lines and cubically-exceptional points as 
an exceptional line. 
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[a-M-l±2A]. 




a± = 



2M+2 
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Figure 2) An enlarged view of the M = 3 phase diagram, showing the 
(a + ,a_) coordinates. The boxes indicate the locations of the quadratic 
and cubic exceptional points, at (a+,a_) = (0,1/2) and (1/4,1), which 
are discussed later in the main text. 

The exactly-zero energy levels can be used to control the pairing-off of eigenvalues and the 
associated formation of exceptional points [2]. On the 'supersymmetric' lines a± = there is 
always at least one zero-energy eigenvalue, for any value of M. The points where this eigenvalue 
becomes degenerate with a second one can be found by looking for zero eigenvalues of the 
supersymmetric partner potential, the partner for (a+,0) being (a+ — , —1) and that for 
(0, a_) being (—1, a_ — ffej)- This idea was used in [2] to show the existence of quadratic 
exceptional points for 

(a+ , a_) = (0, m - j^) and (a+ , a_) = (m - jfa , 0) (2.2) 

where m E N = {1, 2, . . . }. At M = 3 these are the points on figure [2] where the cusped curve 
touches the lines a± = 0. For M = 3, a similar argument can be applied on the other lines 
a± = n £ N on which there is an exact zero-energy level, using a higher-order supersymmetry 
to eliminate this level together with 2n others [1,2]. This establishes the existence of quadratic 
exceptional points at 

(a+ , a_) = (n , m — i) and (a + , a J) = (m — | , n) , m£N,nG Z + (M = 3) . (2.3) 

These are the points on figure [3] where cusped curves touch the other lines a± G Z + . In the 
next section we will generalise these results to other values of M. 

2.2 Locating exceptional points using self-orthogonality 

Our alternative argument starts from the idea, discussed in, for example, [16], that at an excep- 
tional point at least one state will be self-orthogonal, in the sense that its inner product with 
itself under a suitable symmetric inner product must vanish. For the present paper we take this 
inner product to be 

{f\g) = Jj{x)g{x)dx (2.4) 

where the contour C is as in section 1. This inner product is bilinear rather than sesquilinear, 
and - at least in cases where the contour C is the real axis - it is sometimes referred to as 
the c-product [18]. Correspondingly we will refer to yJJ\f) as the c-norm of /; note that 
there is no need for this to be a real number. The c-product is well-defined for any pair of 
functions which decay exponentially as \x\ — > oo along C, and Ti is symmetric with respect to 
it: (f\Hg) = (Hf\g). 
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At an exceptional point, associated with some eigenvalue E, the Hamiltonian acquires a 
Jordan block form, and a so-called Jordan chain {^>(°) . . .ip^^ 1 ^} can be defined which spans 
the subspace of the k merging levels, such that 

(H-E)^=^- x \ j = 0...fc-l, V M) = 0. (2.5) 

Then (V> (0) |V> {0) ) = (^\(H-E)^) = {(Ti-E)^^) = 0, and so the state V (0) is indeed 
self-orthogonal with respect to the c-product. Conversely, suppose that some eigenstate ip, with 
eigenvalue E, has vanishing c-norm, so that (V'lVO = 0- We would like to show that this implies 
that our system is lying at an exceptional point, and to this end we recall a useful result, 
previously exploited in this context by Trinh [19]. Suppose that E n is an eigenvalue, so that y_i 
and yi are proportional to each other and T(E n , a, A) = 0. In fact, from (|1.8p . for such an E we 
have y-i(x, E,a, \) = —yi(x,E,a,X) . Writing ip = y_i = —yi, the relevant result, converted 
to the normalisations used in this paper, is 

m)\ En =T'{E n ) (2.6) 

where the prime denotes differentiation with respect to E, and the dependence of T on a and A 
has been left implicit. Now suppose that, as a function of some combination of a and A, 
has an isolated zero. Then at this point, T(E) as a function of E has a multiple zero which it 
does not possess in the neighbourhood of this point. The zero of (V'lV') rmist therefore mark a 
point where two or more eigenvalues of the eigenproblem have collided. Given the impossibility 
of genuine degeneracies in this problem, this must be an exceptional point, as claimed. 

These results are useful in the present context because along the lines a+ = n and a_ = n, 
n £ Z + , one eigenfunction can be found exactly, namely that with eigenvalue E = 0. Consider 
the line a_ = n, along which 

a-2A = (2n+l)(M+l). (2.7) 

(Corresponding results for the line a+ = n can be obtained by negating A throughout in the 
following.) Then the zero-energy eigenfunction ip = y_i = —yi, normalised in line with (|1.4[) 
and (jl.6p . is 

_ 1 n! (M+l)° , +> -2 W M + - x (a)W/( M +1) 

y 1 ^ 2™ y ' V M+l J y ' 

where L\P (t) is the n th generalised Laguerre polynomial. (To check that ip has been normalised 
with the correct asymptotic, note the relation (|2.T|) between A and a and the fact that the 
highest term of Ln\t) is ^3 t n .) 

To evaluate (ip\ip), we distort the contour C to the union of rays —7-1 + 71, where 

7±i = {x = le^ M+1 \ t e [0, 00)}, (2.9) 



and then use the integral (|C.1|) . analytically continuing in A if necessary to ensure convergence. 
The final result is 

(V#) = £ (M±iy—Tm _ t__ Q n (\) (2.10) 



r(i-M+i, 



where Q n (A) is a polynomial of degree n in A, which can be expressed in terms of the hyperge- 
ometric function 3^2 and Pochhammer symbols (x)^ = x(x+l) . . . (x+k— l)fe as 

Qn(A) = (l-j^) n (l+jj^j) n3 F 2 (-n, ||^f, jf+j', ~ n +TT+T> ^) 

E("l) fc (^) (1-1^1-^(1^)^1+^+^ • (2.11) 



fc=0 
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The zeros of (|2.10j) locate all those exceptional points on the line a_ = n which involve the 
merging of levels at the eigenvalue E = 0. For M = 3, we will argue in the next subsection 
that this captures all exceptional points on this line, with the zeros of Q n (A) being the cubic 
exceptional points associated with cusps on the phase diagram. For other values of M, as will 
be described in more detail in section [H the cubic exceptional points move away from the lines 
a± = n, to be replaced on these lines by pairs of quadratic exceptional points, only one of each 
pair being at E = and corresponding to a zero of Q n (A). 

By contrast, the infinitely-many zeros of the factor l/r(l— f^rj ) in (|2.10j) always correspond 
to quadratic exceptional points. These zeros are at 

2A = (M+l)m-2, m G N (2.12) 

and using (|2.7p they imply the existence of exceptional points at 

(a+,a_) = (n + m— j^px,n), m G N , n G Z + . (2.13) 

This result matches and extends the previously-known cases: for n = 0, it yields the points f|2.2j) . 
found in [2] using ideas based on super symmetry, while for M = 3 the result (|2.3j) is reproduced. 



2.3 Locating exceptional points using quasi-exact solvability 

Self-orthogonality yields important information about the phase diagram at general M, but it 
fails to identify the degrees of exceptional points, and it only sees exceptional points which have 
eigenvalue zero. In this subsection we describe a complementary tactic, special to M = 3, which 
avoids these problems by exploiting the fact that for M = 3 the model is quasi-exactly solvable 
(QES) on the lines a± G N. This will allow us to prove some general statements about the 
spectrum of the model on these lines. A key part of the argument, established in [2], is that any 
complex levels on the lines a± G N must lie in the QES sector of the model. 

For the rest of this section and all of the next we therefore restrict to M = 3, and, to 
minimise the proliferation of factors of i, we replace x by z = ix and set $(z) = tp(z/i). The 
quantisation contour is also rotated by 90°, and the eigenproblem (II. lj) becomes 



d a 2 A 
+ z° + az 2 + 1 



2 _ 1., 

'&(z) = -E$(z), $(z) G L 2 (iC). (2.14) 



dz 2 z 2 

A choice for the contour iC which avoids all singularities in the wavefunctions is given in equa- 
tion (|3.3p below; alternatively a rotated version of ([2.9]) can be used, with suitable analytic 
continuations whenever divergent integrals are encountered. 

If boundary conditions had been imposed at z = and z = +oo, the problem (|2.14j) would 
have been quasi-exactly solvable whenever a and A were related by a = —(4 J + 2 A) for some 
positive integer J, with J energy levels exactly computable [20]. Bender and Dunne [21] found 
an elegant method to find the corresponding wavefunctions, square integrable along the positive 
real axis. We are instead interested in solutions defined along the contour iC, but with minor 
modifications the approach of [21] can still be used3. We set J = q/4 — A/2 and look for solutions 
of the form 

= e^z x+1 2 a n (X)p n (E, A, J) z 2n (2.15) 

n=0 

where 

- (A)= (-i)" ^r(nlA + l) - (216) 
The function 3>(z) will solve (|2.14]) if the coefficients p n (E, A, J) satisfy the recursion relation 

p n = -E Pn - 1 + 16(J-n + l)(n-l)(n-l + \)p n -2, n>\. (2.17) 



^For earlier discussions of quasi-exactly solvable 'PT-symmetric sextic potentials, see [1,2,22]. 
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Setting po = 1 fixes the normalisation, and then p\ = —E follows from (|2.17p at n = 1. If J 
is a positive integer, then the second term on the RHS of (|2.17p vanishes when n = J + 1, and 
so p,j+i is proportional to pj, as are all subsequent coefficients p m> j + i. At a zero of pj the 
series therefore terminates. Owing to the sign of the argument of the exponential prefactor in 
(|2.15p (opposite to that in [21]), the corresponding &(z) will automatically satisfy the revised 
boundary conditions. We define the J th Bender-Dunne polynomial for this problem to be 

Pj(E,X)=pj(E,X,J) . (2.18) 

This is a polynomial of degree J in E, and degree J — 1 in A. By the above reasoning, its zeros in 
E give the J quasi-exactly solvable (QES) levels that the model possesses on the line a = 4J+2A. 
Since boundary conditions are not imposed at the origin, replacing A by —A throughout also leads 
to an acceptable solution, and so for each J £ N there are two lines of quasi-exact solvability in 
the (2A, a) plane: a = 4 J + 2A and a = 4 J — 2A. In the (a+,a_) coordinates these lines are 
(a+,a_) = (|( J— 1+A), t^( J— 1)) and (^( J— 1), ^( J— 1— A)) respectively. Figure[6l below, shows 
the QES lines on the (2A, a) plane. 

These lines are very useful in mapping the exceptional points on the whole (2A, a) plane. 
The reasoning is best explained via a sequence of lemmas, which may be of independent interest. 

Lemma 1: The Bender-Dunne polynomials satisfy the 'reflection symmetry' 

Pj(E,X) = (-i) J Pj(iE,-J-X) (2.19) 

Proof: Introduce a set of polynomials defined by r n (E, A, J) = (—i) n p n (iE, — J— A, J). Direct 
substitution into (|2.17p shows that the r n satisfy the recursion 

r n = -Er n - 1 + 16(J-n + l)(n-l)(J-n + l + X)r n - 2 , n>l (2.20) 

with initial conditions r_i = 0, r$ = 1. The claimed symmetry is equivalent to rj(E, X, J) = 
Pj(E, A, J). Now consider a more general recursion 

u n = -Eu n -i + 6„_iu„_ 2 , n > 1 , u-i = , u = 1 ( 2 -21) 

with some set of coefficients {b n }. It is straightforward to verify that the general solution is 

[n/2] 



[n/2j / \ 

E b n b t2 ...h k \E n - 2k . (2.22) 

fc=0 \{0<i 1 <--<i k <n} ) 



In particular, pj is given by (|2.22p with n = J and hi = 16( J— l)(i— 1+A), and rj by 
(|2.22p with n = J and bi = 16(J— 1)(J— A). Thus the two differ by the substitution 
fri - * bj-i, and since (|2.22p for n = J is itself symmetrical under this mapping, the lemma is 
proved. 

Lemma 2: If A < 1 — J, then all zeros of Pj(E, A) are real and distinct. 

Proof: The given values of A correspond to the points on the QES line a = 4J + 2A which 
lie in the region a < 4 + 2|A|. The reality result (|1.9j) then implies that the spectrum of the 
eigenproblem (|2.14p . which includes the QES sector described by the zeros of Pj(E,X), is real. 
These zeros must therefore all be real. To show that the zeros are simple, we use the fifth 
spectral equivalence from [1] to map our problem on to one which, for the given range of A, is 
hermitian. Converted into the current coordinates, this equivalence states that the spectrum 
of (|2.14p is the same as that of the following radial problem for functions defined on M + and 
decaying at x — > +oo: 



r/2 \'2 _ I 



dx 2 x 2 



<j>(x) = E (j>(x) , <f>(x)\ x ^ ~ x x ' +1 ' 2 , (2.23) 
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where A' and a! are related to A and a by 




(2.24) 



This equivalence is illustrated in figure[Sl the dual problem on the right is hermitian for A'+l/2 > 
— 1/2, which translates into 2A + a < 4 on the left. The QES line under discussion is a = 4J+2A; 
it maps onto the diagonally-oriented QES line a 1 = — 4J — 2A' on the right-hand diagram, and 
is in the hermitian region of that plane for A < 1 — J. Since all eigenvalues of the hermitian 
spectral problem are distinct, so must be the zeros in E of Pj(E,X), for all A < 1 — J. (It is 
worth remarking that the standard proof of the simplicity of the eigenvalues for the hermitian 
problem uses essentially the same steps as lead to (12. 6p , together with the fact that for hermitian 
problems the eigenfunctions can be taken entirely real, so that the LHS of (12. 6p never vanishes.) 



K (0.4J) 



'HJ.O) 



(4J.0) 



2X 



2X\ 



(-2J.6J) + 



(-2J,-2J)X 



a / 



a' 



2X 



\(2J,-6J) 



Figure [5] The QES lines a = 4 J + 2A and a = 4 J — 2A in the (2A, a) plane, and their images 
in the (2A', a') plane under the mapping (I2.24j) . On the right, the images of the 2 A and a axes 
are also shown. The left-hand diagram corresponds to lateral boundary conditions, while those 
for the diagram on the right are radial, and hermitian when 2 A' > — 2. 

Aside: To make the proof of lemma 2 self-contained, it is possible to show directly that Pj(E, A), 
as defined for the non-hermitian problem (I2.14h . is equal to the standard Bender-Dunne poly- 
nomial for the Hermitian problem (|2.23p . QES levels for (|2.23[) occur when J' = — (a' + 2A')/4 
is a positive integer, or in other words when a' = —4 J' — 2A'. Defining p' n (E, A', J') to satisfy 
the recursion 

p' n = Ep' n __ x + 16(n — l)(n — J' — l)(n + A' - \)p' n _ 2 , (2.25) 

with p'q = 1, the QES levels are given by the zeros of the polynomial Pj,(E, A') = p'ji[E, A', J') 
(see [21] for details). If A' and a' are given in terms of A and a by (|2.24p . then (|2.25j) becomes 

Pn = Ep' n ^ + 16(n - l)(n - J - l)(n - J - A - l)p' n _ 2 (2.26) 

where J = (a — 2A)/4 = J'. Since this recursion is, up to a swap E — * —E, the same as (I2.20p . 
it follows from the proof of lemma 1 above that P'j = Pj, as claimed. 

Lemma 3: If A > —1, then all zeros of Pj(E,X) are purely imaginary (or zero) and distinct. 
Proof: This follows from lemmas 1 and 2. 
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Lemma 4: On the infinite segment A > — 1 of the QES line a = 4 J + 2 A, the eigenproblem 
f|2.14j) has exactly J distinct non-real (in fact purely imaginary) eigenvalues for J even, and J — 1 
non-real (and purely imaginary) eigenvalues and 1 zero eigenvalue for J odd. By the A — > —A 
symmetry of the problem, the same holds for the A < 1 segment of the a = 4 J — 2 A QES line. 
Proof: First recall from [2] that on QES lines, the non-QES sector of the spectrum is entirely 
real. The result then follows on combining lemma 3 with the fact that the degree J polynomial 
Pj(E, A) is a function of E 2 for J even, and E times a function of E 2 for J odd. 

The results so far show that on the QES lines a = 4 J + 2A the spectrum is entirely real for 
A < — J + 1, and has exactly 2 [J/2] complex eigenvalues for A > —1, where [x] denotes the 
largest integer less than or equal to x. This is illustrated in figure [6j 




V 

-15 -10 -5 5 10 

IX 

FigureO The lines a = 4J±2A, J € N, of quasi- exact solvability on the (2A, a) plane. 
The arrows indicate the (open) segments of the lines a = 4J+2A on which the precise 
numbers of non-real eigenvalues are determined by lemmas 2 and 4. Note that these 
results are consistent with the locations of the numerically-obtained curved cusped 
lines (blue online) , across each of which the number of non-real eigenvalues increases 
by two. For J ^ 1, only those parts of the QES lines which lie in the zone of possible 
unreality a > 4 + 2|A| are shown. The lines with J odd coincide with the protected 
zero-energy level lines shown in figures [2l [3] and [4j 

The situation in the remaining intervals — J + 1 < A < —1 is clarified by lemmas 5 and 6. 

Lemma 5: At the points A = — J + n, n = 2 . . . J— 1 on the QES line a = 4J + 2A, the problem 
has exactly 2[n/2] distinct non-real eigenvalues. 

Proof: In addition to being on the line a = 4J + 2A, the given points lie on the lines a = An — 2 A 
in a region where lemma 4 applies. 

Lemma 6: On each segment — J + 2m — 1 < A< — J + 2m of the QES line a = 4 J + 2A, where 
J > 2 and m = 1, 2 . . . [(J+l)/2] — 1, there is at least one point where the eigenproblem has an 
exceptional point with eigenvalue zero. The same is true of the segments J— 2m < A < J— 2m+l, 
m = 1, 2 . . . [(J+l)/2] - 1 of the QES line a = 4 J - 2A. 

Proof: By lemma 5, when A = — J + 2m — 1 the number of non-real eigenvalues is 2m — 2, while 
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when A = — J + 2m it is 2m. The number of non-real eigenvalues thus changes by two as A 
moves from — J + 2m — 1 to — J + 2m. By the VT symmetry of the problem, non-real eigenvalues 
always occur in complex-conjugate pairs; combining this with the E — > —E symmetry of the QES 
sector, any non-real eigenvalues created or destroyed away from E = must appear in quartets. 
So to change the number of non-real eigenvalues by two, at least one pair must be created or 
destroyed at zero, and hence there must be at least one exceptional point with eigenvalue zero 
in each interval — J + 2m — 1 < X < — J + 2m. The final statement of the lemma then follows 
from the A — > — A symmetry of the problem. 

Lemma 7: For M = 3, the zeros of the polynomials Q n (X), defined in equation (|2.1ip above, 
are all real and simple, with one in each interval A E [2m— 1, 2m], m = 1 . . . n. 
Proof: Specialising the discussion of subsection 12.21 to M = 3, Q n {\) has a real zero at every 
point on the line a = 4(2n+l) + 2A, A 6 1 where the eigenproblem has an exceptional point 
with eigenvalue zero. Combining this with lemma 6 taken at J = 2n + 1, Q n (A) has at least one 
real zero in each interval [2m— 1, 2m], m = 1 . . . n ; but since Q n (X) is a polynomial of degree n 
this must exhaust all of its zeros, which must therefore also all occur singly. 

Lemmas 6 and 7 show that on the lines a = 4 J ± 2 A with J = 2n + 1, there are n odd-order 
exceptional points with eigenvalue zero. If all of these have the lowest possible degree, that is 
three, then at each a pair of complex eigenvalues is created, and would account precisely for 
the total number of complex QES levels which must appear as the QES line is traversed. In 
principle, one could imagine more complicated scenarios where further pairs of complex levels 
appear at these exceptional points and then annihilate with each other later, but our numerical 
results show no evidence of such behaviour and we shall proceed on the assumption that it does 
not occur. If we further assume that the triply-exceptional points are isolated in the sense that 
there are no other triply-exceptional points in their immediate neighbourhoods in the (2A, a) 
plane, then these points must be occurring where two lines of double degeneracy meet at a cusp, 
as illustrated in figure [7] below. Thus, subject to the two assumptions just mentioned, we have 
proved that for M = 3 the cusps in the exceptional lines do indeed lie on the lines of protected 
zero-energy levels, as conjectured earlier. 




Figure [7J The behaviour of the energy level surface E(a + , a J) 
in the vicinity of a cubic exceptional point. 
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The above results capture the [J/2] transitions to complex levels that occur as the QES 
line a = 4J + 2A is traversed, at each of which a pair of complex levels is created. However 
this does not necessarily exhaust all of the exceptional points on the corresponding QES line 
- indeed, for J odd (I2.12p shows that there are infinitely- many more exceptional points on the 
line a = 4J + 2A, at (2A,a) = (4m— 2, 4J+4m— 2), m £ N. Examination of figure [6] reveals that 
at each of these points an exceptional line touches, but does not cross, the QES line, so that the 
exceptional point does not cause the creation of further complex levels while motion is restricted 
to the QES line. Once the QES line is left, the J QES levels start to mix with the non-QES 
sector, and further complex levels can be formed. (In fact, if the point (4m— 2, 4J+4m— 2) on 
the QES line a = 4 J + 2A is left in a direction perpendicular to that line, a pair of complex 
levels is created in a different QES sector, that for the QES line a = 4(J+2m— 1) — 2A.) This 
general picture, and also the claimed isolation of the triply-exceptional points as illustrated in 
figure will be supported by some perturbative calculations in the next section. 

One would also like to be able to rule out the more exotic scenarios for the behaviour of 
levels in the QES sector, mentioned in the discussion following lemma 7. We do not have a 
rigorous argument for this, but we do have extensive numerical, and some analytical, evidence 
for the following conjecture which, if true, would eliminate such possibilities: 

Conjecture: For all A, the squared zeros (in E) of the Bender-Dunne polynomials Pj(E,X) 
are real, and, apart from the zero at E = when J is odd, they are monotonically-decreasing 
functions of A. 

Immediate consequences of a proof of the conjecture would be that QES levels, once complex, 
remain so, and that the only way that QES levels can become complex is via E = 0. In turn this 
would prove that the zeros of the polynomials Q n {X) do indeed correspond to triply-degenerate 
points in the spectrum, since one can easily rule out the presence of zero eigenvalues in the 
non-QES sector at the relevant points. 

The conjecture is similar in spirit to the Feynman-Hellman theorem, but the eigenproblem 
here is not hermitian, and this invalidates any variant of the standard proof. Note also that the 
conjecture is certainly false for the non-QES part of the spectrum, the (un-squared) levels of 
which can pass through zero while remaining real. As a sample of our numerical checks, figures 
[U and [9] show the squared QES levels for J = 20 and J = 21. Apart from the E 2 = line 
on figure all the curves are monotonic; given this, the fact that they all pass through zero 
between A = 1 — J and A = — 1 follows from lemmas 2 and 3 above. 




— i 1 1 1 1 ■ r 

-30 -20 -10 



X 

Figure [5] Squared QES levels for J = 20. The dotted vertical lines are at A = 1 — J and 
A = — 1 ; all transitions from real to imaginary eigenvalues occur for 1 — J<A<— 1. 
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Figure O As figure [51 but for J = 21. Again, all transitions to imaginary eigenvalues 
occur for 1 — J < A < — f. 

The curves shown on figures [8] and [9] appear to asymptote to linear functions of A as A — > ±oo. 
This turns out to be the will be shown below, where the slopes of these functions 

will also be found exactly. Figures [10] and [11] illustrate these results, and further support the 
monotonicity conjecture, by plotting the derivatives of the squared QES levels, again for J = 20 
and J = 21. 




-30 -20 -10 



I 

Figure[l0l Derivatives of squared QES levels for J = 20. The straight horizontal lines (red 
online) show the predicted asymptotic values for these derivatives, which are everywhere 
negative. The symmetry of the plot about A = J/2 follows from the reflection symmetry 
proved in lemma 1. 
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Figure [TT1 As figure [TU1 but for J = 21. All but one of the plotted functions are negative, 
the exception corresponding to the level at E = 0. 

To treat the large- |A| behaviour of the QES levels analytically, a first approach is to return 
to the Bender-Dunne recursion relation (|2.17p . In the limit |A| S> J, and keeping n < J, this 
simplifies to 

Pn = -£pn-i + 16A(J -n + l)(n-l)p n _ 2 , Po = l- (2.27) 

Solving for low-lying values of J, a remarkable simplification occurs precisely at n = J, where 
the asymptotic Bender-Dunne polynomials pj- symp are found. For J even, 

J/2 

jP asymp ( ^ ) = JJ + 16(2fc-l) 2 A) , (2.28) 

fe=l 

while for J odd, 

(J-l)/2 

Pp mp (E) = -E Yl (E 2 + 16(2/t) 2 A) . (2.29) 
fc=i 

Hence the squared QES levels are indeed linear functions of A in this limit, with slopes which 
are negative, and proportional to the squares of odd or even integers. Rather than prove these 
formulae directly from the asymptotic Bender-Dunne recursion relation, we will use the same 
spectral equivalence as employed in the proof of lemma 2 above. The polynomial Pj(E,X) 
encodes in its zeros the QES levels along the line a = 4J + 2A; by the reflection symmetry (12.191) 
it will suffice to consider the asymptotic behaviour of the levels in just one direction, which we 
choose to be A — > — oo. Then, instead of applying the mapping (12.240 immediately, we precede it 
by the trivial symmetry (2A, a) — ► (— 2A, a) of the "PT-symmetric problems, which flips between 
the two QES lines shown on the left-hand diagram of figure 03 The line (2A, 4 J + 2A), A £ M, is 
now mapped to (— 2J, — 4A) on the (2A', a') plane, the vertical line on the right-hand diagram of 
figureEJ The limit A — ► — oo is thus mapped to a' — > +oo, 2A' = — 2J in the spectrally-equivalent 
lateral problem (|2.23p . For the QES sector, we are only interested in the J lowest-lying levels, 
where J remains fixed as a' = — 4A — * +oo. In this limit the quadratic term in (I2.23h comes to 
dominate and the problem reduces to the (scaled) radial simple harmonic oscillator 

l 

4>{x) = E4>{x), <A(x)U_ ~x- J+1 / 2 . (2.30) 



dx 2 



\\x A 
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Since J > 1, the boundary conditions at the origin are irregular. Nevertheless, the problem can 
be solved exactly for any value of J, the first J levels being 

E = \/— 4A(— 2J + 4n - 2) , n = 1,2, ... J. (2.31) 

(Strictly speaking, a resonance when J is an integer means it is safest to shift J slightly away 
from integer values for the calculation, but the final result is unaffected.) It is straightforward to 
see that this confirms equations (|2.28[) and (|2,29p above, as desired. (The overall normalisation 
of p^ sympt (£;) can be fixed by considering the coefficient of E J .) Note that these results imply 
the truth of the monotonicity conjecture in the limits |A| — > oo and thus give some supporting 
evidence for its general validity. Unfortunately, the regions |A| — > oo are not of interest from the 
point of view of reality properties, since they are already covered by lemmas 2 and 3 above, and 
so a full proof of monotonicity of the squared eigenvalues would still be worthwhile. 

Finally, an alternative way to locate the cusps corresponding to the collision of levels at E = 
is to examine the odd Bender-Dunne polynomials fbm+l (E, A) directly. These factorise as E 
times a polynomial in E 2 , and there will be a multiply-degenerate zero-energy level whenever 
this polynomial vanishes at E = 0, or equivalently whenever 

Jip 2m+1 (E,\)\ E=0 = Q. (2.32) 

For fixed m and J = 2m + 1, consider the sequence of polynomials p2 n +i(E, A, J), n = 0. . . m, 
and define 

q n (X,m) = ^^ —p 2 n+i(E,X,2m+l)\ E=0 . (2.33) 

A consideration of the Bender-Dunne recurrence (|2.17p and its derivative at E = shows that 
q n (X,m) satisfies the first order recurrence 

q n = (m-n + ±)(n + ±\)q n -i + (™) f[(k - §)(* - \ + \\) , (2.34) 

k=i 

with initial condition qo = 1. Anticipating the final result in our notation, we set 

Qm(A) = q m {X,m) (2.35) 

so that Q m (X) is a polynomial in A of degree m, and its zeros are the points identified by (|2.32p . 
For example: 

Qi(A) = ^(3 + 2A) (2.36) 

Q a (A) = -3-(41 + 40A + 8A 2 ) (2.37) 
lb 

QsW = 7^(7 + 2A)(63 + 56A + 8A 2 ). (2.38) 
64 

It turns out that the general solution to (|2.34|) can be expressed using the hyper geometric 

function 3F2 and Pochhammer symbols (x)k = x(x+l) . . . (x+k— 1), {x)k = (— l) fc (l— k— x)k ■ 
For n = m this solution is 



QmiX) = = {l+\X) m {\) mZ F 2 {-m, \,\+\X; l+±A,i-m; 1) 

EC" 1 )* (Tj (l-^)-(|+|A) fc (l+iA + fc) m _ fc . (2.39) 



This is the M = 3 case of the general formula (|2.1ip . here derived by a completely different route. 
Since terms being summed in (|2.39[) are invariant up to a factor of (— l) m under the simultaneous 
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4 111 8 ' 4 111 8 



9 9 i \/l70+772li 9 i A/l70-7y / 2jf 
4> 4 =t g ! 4 =t g 



Table 1: Location of some of the cusps in the («+, a_)-plane for M — 3. 



exchange k — ► m—k, A — > — 2m— 1— A, the M = 3 polynomial Q m (A) is invariant up to a sign 
under A — > —2m—l—\ and its zeros are symmetrically distributed about A = — m — 1/2. This 
reflects the more general result (|2.19p . For m odd, this means that Q m {—m — 1/2) = 0. 

Returning to the a± coordinates, if J = 2m+l and m G N then the QES line a = 4J + 2A 
corresponds to (a+, a_) = (m+A/2, m), and the m cusps on this line occur at the zeros of Q m (A), 
while for a = 4J — 2 A the cusps lie on the line (aq_, a_) = (m, m— A/2) with A a zero of Q m {— A). 
The first few cases from this second set are given in table HI to within our numerical accuracy, 
they match the cusp positions shown in figure Notice that the symmetrical distribution of 
the zeros of Q m (X) mentioned at the end of the last paragraph implies a relationship between 
the locations of pairs of a priori unrelated cusps and gives a simple formula for the remaining 
'unpaired' cusps: (a+,a_) = (m,m/2— 1/4) and (a+,a_) = (m/2— 1/4, m) for all odd m 6 N. 



3 Jordan blocks for M = 3 

3.1 The Jordan block at a quadratic exceptional point 

We now investigate the exceptional points and their neighbourhoods in more detail, beginning 
with a quadratic example. The first step is to find the Hamiltonian Hq at the exceptional point, 
restricted to the two-dimensional space of states which merge at that point. This will have a 
Jordan block form. We then perturb about this point by writing the full Hamiltonian, H, as 
H = Hq + V, and expand H using the wavefunctions of Hq as a basis. Thus we will need to 
calculate 

Hmn = {(t>m\H\(j) n ) = (4> m \H \(j) n ) + {^ m \V\(j> n ) , (3.1) 

where {cp m ,4>n} 1S a basis for the Jordan block form of Hq, and {4>m,4> n } is an appropriate 
dual basis, so that the functions together form a part of a biorthogonal system, as discussed 
in [23] for generic cases and [16] in the presence of exceptional points. In the current setting, 
wavefunctions decay as \z\ — ► oo along iC and a suitable pairing between functions g(z) and 
'dual functions' [23] f(z) is a rotated version of the c-product ([22 



</l<?) = / f(z)g(z)dz . (3.2) 
he 

Here and below a convenient choice for iC, beginning and ending in the (rotated) Stokes sectors 
iS-i and iS%, will be 

iC = -7_i + 70 + 71 (3-3) 

where j±i = {te ±nt ^,t G [s, oo)}, 70 = {ee lt ,t G [— rr/4, 7r/4]}, and the small positive number 
e ensures that any singularities at z = are avoided. For later use we note the following basic 
integral along the contour iC, which holds for arbitrary o£l: 

o(o-3)/4_,- 

fS' 2 dz= 2 " ™ . (3.4) 
iC r(i(3-a) 1 7 
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This is easily checked, either by analytic continuation from a > — 1 or using the real integral 
J°° t a e~ l l 2 dt = 2^ a ~ 7 ^ 4 r(|(a+l), ^e 4 ), where T(a,z) is the incomplete gamma function. 

As our example we take the exceptional point at (2A, a) = (2, 6), (a+, a_) = (1/2, 0), shown 
with a box on figured! This lies on two lines of quasi-exact solvability, (a = 4 J + 2A)|j=i and 
(a = 4 J — 2A)|j=2 , and we shall focus on the second of these. Along this line a = 8 — 2 A and 
the eigenvalue problem is 

A 2 - - \ 

- 2X)z 2 + + E U = . (3.5) 

Setting A = 1 — 2e, this corresponds to (a+, a_) = (1/2, e) and the exceptional point, where two 
levels merge and the Hamiltonian can be written in a Jordan block form, occurs at e = 0. The 
recursion relation (|2.17p for p n {E, — A, J) becomes 

Vn = -Evn-\ + 16(3 - n)(n - l)(n + 2e - 2)p n „ 2 (3.6) 

and, as expected, the second term on the RHS vanishes when n = J + 1 = 3. The energy 
eigenvalues of the two QES levels are given by the roots of P2(E,— A, 2): E?o,i = ±4z\/2e = 
±4\/ — 2e. The corresponding (unnormalized) eigenvectors are, from (|2.15p . 




% = ^VS(lT^)ef (3.7) 

where ^0,1 = ^p^^ cq ' ^'^ \E 1 ■ At e = these two eigenvectors coincide, and to see the 
Jordan form of the Hamiltonian we proceed as in appendix |A] and construct 



^0 = *0| e=0 

= aV2z- 1 / 2 e z4 / 4 (3.8) 



de 



e=0 

where a is a constant. The Hamiltonian at e = is 

Ho = -^ + z 6 + 6z 2 + ^ (3.9) 

and requiring 0o and 0i to satisfy the 'Jordan chain' relations Hq 0o = and Hq 0i = 0o fixes 
a = — 1 an d shows that the Hamiltonian has the desired Jordan block form. However, this 
basis is not unique: replacing {0o,0i} by {jU0o,A*0i + ^00 } preserves the Jordan chain for any 
constants /J, and v. This freedom can be used to make a convenient choice for our biorthogonal 
system. Dropping a common factor of —i, the general Jordan basis is 

= ^ 3 /V 4 / 4 (3.10) 



The integrals f iC c, 


W 


) n & with m,n £ {0, 1} can be evaluated using 1 


3.4p and are 


/ 0000 




= /x 2 / zV 4 / 2 dz = o 




/ 0001 
JiC 




= / (l^ z + ^) e ^dz= l -^V2^ 
JiC ° 




/ 0101 




= [ ^ z -l-y uz + u 2 z ^ e * 4 /2 dz= it 
JiC 4 


( g7T// + a/27TZ^ 


Therefore, if 









(3.12) 
(3.13) 
(3.14) 



4V2 V '7T 
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then J iC 4>\(f)\ dz = and f iC 4>o4>i dz = 1, allowing us to take the dual basis to be <fio = <p\ and 

4>i = <p . 

The spectrum in the neighbourhood of Hq can now be investigated. For this a two-parameter 
family of perturbations is required, and we take one of these parameters to be e, and the other, 
77, to be orthogonal to this in the (a+,a_) coordinates so that (a+,a_) = (1/2 + 77, e) and 

H = Hq + 4(r/ + e)z 2 + 4(1 + 77 - e)(r/ - e)z" 2 . (3.16) 

The matrix elements of interest are, in an obvious notation, 

<0o,ik 2 |0o,i) 

(0o,i \z~ 2 1 <A),l) 
and so the truncated Hamiltonian, to leading order, is 

*W « ( 16^326 2^77 ) • (3 - 19) 

The approximate energy levels are thus the roots of the characteristic polynomial of this matrix: 

E = 2^77 ±4vV - 2e. (3.20) 

Two special cases deserve comment. For 77 = the QES levels E = ±4\/— 2e are recovered, as 
expected given that these QES levels were used to set up the approximation scheme in the first 
place. More interesting is the fact that for e = the energy levels remain real as the exceptional 
point is traversed. (The same phenomenon was remarked in a finite-dimensional setting in [15] .) 
This corresponds to the fact that the approximation correctly predicts the direction of the line 
of exceptional points away from the QES point (a+,a_) = (1/2,0). However one should be 
wary of trusting the approximation any further - one might expect that the curvature of the 
line of exceptional points could be recovered from the line of points where the discriminant of 
the characteristic polynomial of f|3.19j) vanishes, which is T7 2 — 2e = , or a_ = | (| — a+) . 
However, a fit to the numerical eigenvalues of the full equation shows that the shape of the 
curve of exceptional points near to (1/2,0) is rather given by a_ ~ k (i - a+) with k « 0.78. 
Given that this curvature is controlled by sub-leading effects, this failure should not be too 
surprising, but it does highlight the delicacy of perturbation theory about exceptional points. 
A more systematic investigation of this issue would be valuable, but for now we will pass on to 
an examination of a typical cubic exceptional point. 



2tt 
1 



4 




(3.17) 
(3.18) 



3.2 The Jordan block at a cubic exceptional point 

From table [H the first cubic exceptional points occur at (a+,a_) = (1,1/4) and (1/4,1), on 
the J = 3 QES lines a = 12 — 2A and a = 12 + 2A. We focus on the line a = 12 + 2A and set 
A = 2e — 3/2 so that (a + ,a_) = (e + 1/4, 1) and the exceptional point occurs at e = 0. The 
eigenvalue problem in terms of e is 

~ + z G + (46 + 9)z 2 + (2£ ~ ~ 2) + * = (3.21) 

and the recursion relation for p n (E,2e — 3/2, 3) is 

p n = -E Pn ^ + 16(4 - n)(n - l)(n + 2e - 5/2)p n _ 2 . (3.22) 
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The roots of give the energy eigenvalues of the three QES levels: Eq = 0, E± = ±8y/— 2e. 
The corresponding eigenstates are 

*o = e^z^a (l + (3.23) 

where a is some normalisation to be fixed later. Note that when e = these three QES eigenstates 
merge and we have only one known eigenstate at this point, namely 

V \ e=0 = ae zi /*(- + 2z 3 ). (3.24) 



3.2.1 The Jordan basis 

If we perturb away from the cubic exceptional point along the QES line, the Hamiltonian will 
correspond to a toy model matrix of the form 



L(e) = e/2 1 . (3.25) 




The method that we used to calculate the Jordan basis for the quadratic exceptional point 
in section 13. 1\ explained for n x n Jordan blocks in appendix |Aj does not apply here. This 
is because the matrix considered in appendix |A] would correspond to a perturbation of the 
Hamiltonian along a line perpendicular to the QES line, along which we do not know the 
relevant eigenfunctions analytically. Instead, we will have to find the basis functions for (13.25D 
by solving the Jordan chain constraints directly, to find wavefunctions (f>o, 4>\ and fa that satisfy 

H (j)o = 

Hofa = <P (3.26) 

Ho4>2 = 4>i 

where Hq is the Hamiltonian at the cubic exceptional point: 

Ho = ~& + z6 + 9z2 + | ■ t 3 ' 27 ) 

Note that Hq ^/o| e =o = so we can take 4>q = ^ r o| e =o- Th en solving (|3.26p for <j>\ and <fo, we 
find 

<h = c^ 4 (f + 60 + 2z 3 )) (3.28) 

^( — + - + c( 1 - + 2z* 
16z 2 V* 

with a, b and c constants, arbitrary at this stage. These are the most general solutions to (|3.26p 
that also satisfy the relevant boundary condition, that is square integrability along iC. 
Now that we have a basis, we must find the dual basis 4>q, <j>\ and 02 which satisfies 

c^i^idz = 1, fori = 0,1,2 

, : c 

pi4>jdz = 0, fori ^ j. (3.29) 

iC 
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From [16] we expect the dual basis to be 0o = <j>2, 4>i = 4>i an d 4>2 = <fio and this is supported 
by the fact that f iC 0o0o dz = J iC (j)o4>i dz = and J iC 0i0i dz (x a 2 . Fixing J- c 0202 = and 
J iC 4>\4>2 dz = constrains the coefficients b and c to be 



OTT 



16r(3/4) 2 
q(37r 2 - 8r(3/4) 4 ) 
C ~ 512r(3/4) 4 

Then requiring J iC 0i0i dz = f iC 0o0 2 dz = I fixes a 2 : 



(3.30) 



a 2 = _±_L, (3.3i) 



r(3/4) 

Choosing the root with positive real part for a we have fixed the basis to be 
(1 - i)2 7 / 8 



00 



3 *V4 Q + 2z 3 ) 



8«r 7 +2tt^ (3.32) 

i 6 r(f) 5/2 V* W 

(l-i)2 7 /8 4 /24r(f) 4 + 3vr 2 /3\ 2 ., ., , ., v ., 

-e 2 /4 ^ 16vrr - z + 6ttV - 16r - z 3 



5i2r(f) 9/2 \ - U7 I 4 

with the dual basis 0o = 02 5 01 = 01 and 02 = 0o- 
3.2.2 Matrix elements and the cusp singularity 

We first perturb away from the cusp along the QES line (a+,a_) = (e + 1/4, 1) and write 

H = H + V , (3.33) 

where V = 4e ^ 6£ + Aez 2 is considered as a perturbation of H (|3.27|) . The required matrix 
elements are 

(HV\<f> ) = , (3.34) 

3r(i) 2 

and 

<0l|F|0O> = (02|^|0l) 

8e / 2 4r^ 4 -i2r^ 4 e + ^ 



3F V4 



|) 4 \ V 4 / V 4 , 

64e (3.35) 



to leading order in e. To investigate the shape of the cusp we also need to perturb away from the 
exceptional point in the direction perpendicular to the QES line, i.e. along 77 where a = — 4ry + 9 
and A = 2ry — 3/2, or a + = 1/4 and a_ = 1 — 77. The Hamiltonian is now H = Hq + V + V with 



V = -4 V z 2 + 2 J^lfR , (3.36) 

z z 



and to first order in 77 we find 



WW— 5^" F(|f' <3 - 37) 
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The remaining matrix elements effect the energy levels only at subleading order in rj and so they 
can consistently be ignored. The resulting truncated Hamiltonian is 



/:/„,,, ~ I -64e II. CUS) 

1 1287rr; ' 



7 




r( 3 /4) 

The matrix (|3.38p has the characteristic polynomial X 3 + 128eX + y(3/^}' 2 = 0. Now the curve 
of exceptional points occurs when dX/de — > oo (or equivalently dX/dr] — > oo). Since 

_ -128X 

" 3X2 + 128e ( 3 - 39 ) 

the requirement dX/de — ► oo fixes 

* = (3.40) 

Substituting this into the characteristic polynomial above and restricting to e < gives the 
following relation between r/ and e: 



For e > the relation (|3.40p is not real indicating that there are no exceptional points in this 
region, which matches our numerical results. In terms of the a± notation, a+ = e + 1/4 and 
a_ = 1 — r] so this relation becomes: 



a ^ 1± l^m^ {i/i _ a+)3/2 (3 42) 

which is valid for a_ close to 1 and << a+ < 1/4. 

A comparison between the prediction (|3.42j) for the line of exceptional points in the vicinity 
of the cusp at (a + ,a_) = (1/4,1) and numerical data obtained from a direct solution of the 
eigenvalue problem is shown in figure [T2l The shape of the curve is accurately reproduced. 
In principle the same calculations could be performed for other cusps, though the relevant 
wavefunctions become more complicated. 
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4 Numerical results for 



Having established the existence of quadratic and cubic exceptional points at M = 3, we now 
explore the situation at other values of M. Whitney's theorem for mappings from the plane to the 
plane [24] implies that the fold and cusp singularities (corresponding to the doubly-exceptional 
lines and triply-exceptional cusp points seen at M = 3) are stable, and so the pattern of cusped 
lines must persist, at least while M remains sufficiently close to 3. Recall also that protected 
zero-energy levels lie on the lines a± = n for all values of M. However, away from M = 3 
quasi-exact solvability is lost, and so one of the properties which confined the cusps at M = 3 
to the lines a± = n, namely the symmetry of the set of merging levels under E — > —E, may no 
longer hold. 

Figures [T3l [TH and [TBI show the exceptional lines for M = 2, 1.5 and 1.3. The plots were 
obtained by a direct numerical solution of the second dual form of the eigenvalue problem, as 
described in appendix [Bj 
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Figure fT3l Exceptional lines at M = 2. 

As predicted, the overall pattern remains the same, but the cusps move away from the 
protected zero-energy lines. The points where the outermost cusped line touches the super- 
symmetric zero-energy lines a± = are known exactly, from (|2.2p . As M decreases from 3, 
they move down from the midpoints between = n and = n + 1 along the lines a± = 0, 
as predicted by the formula fl2 . 13H . At the same time, the numerical data shows that the cusps 
move upwards, on the rescaled coordinates of the plots which keep the lines of protected zero- 
energy levels at constant locations. As M — > 1 + the pattern shows signs of simplifying, with the 
cusps heading away towards a = +oo and the regions of unreality shrinking towards the lines 
2A G 2Z. This behaviour will be discussed further in section [5j 
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Figure [TU Exceptional lines at M = 1.5. 
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Figure [TBI Exceptional lines at M = 1.3. 

It is interesting to see the fate of the exceptional points corresponding to the zeros of the 
polynomials Q n (X), which at M = 3 are triply-exceptional cusps. For M/3 the cusps move 
away from the lines a± £ N, and so the zeros of the Q n (^) are n o longer cusps, but are instead 
only doubly exceptional. Furthermore, the presence of an exactly-zero level on the lines a± G N 
forces the smooth parts of the exceptional lines to be tangent to these lines immediately M moves 
away from 3, and this leads to a complicated change in the shape of these curves, illustrated in 
figure [TQ 
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M = 1.5 



M = 3 



M = 6 



Figure [TU The movement of a cusp for M 7^ 3. 

The plots of figure [TU] indicate that for M > 3 the movement of the cusps away from the 
lines a± G N is opposite to that for M < 3, and this can also be seen in figures [T7] and [181 
which show the exceptional lines for M = 10 and M = 30. Again, the locations of the zero- 
energy exceptional points on the lines a± confirm the formula (|2.13p . and there are no hints 
of any further exceptional points beyond those predicted by our general considerations. As the 
cusps move towards the a axis, they start to merge to leave isolated 'islands' of unreality in the 
phase diagram. In the theory of singularities, this merging of two cusps is sometimes called the 
'beaks' transition (see, for example, [25] and references therein). As for M — > 1 + , the structure 
simplifies as M — > oo. 

It turns out that the simplifications near to M = 1 and M = oo can be understood ana- 
lytically, using the fact that the limiting points M = 1 and M = oo are exactly solvable, and 
allow perturbative treatments to be set up in their vicinities. In the next two sections this will 
be developed in detail, starting with the region near to M = 1 where we will see that it leads to 
a novel insight into the original 'Bender-Boettcher' phase transition to infinitely-many complex 
levels, which occurs when M becomes smaller than 1. 
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Figure [T71 Exceptional lines at M = 10. 
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Figure [15] Exceptional lines at M = 30. 

5 Perturbation theory about M = 1 

5.1 Exceptional points via near-degenerate perturbation theory 

In this section we revert to the original formulation of the eigenvalue problem, namely 

H M ip(x) = Eip(x) , if>(x) € L 2 {C) (5.1) 

where 

H M = - (**) 2M ~ a(i X ) M -l + ^1 . (5.2) 

For M = 1 this problem can be solved exactly - it is the "PT-symmetric simple harmonic 
oscillator [11,26], and its spectrum is entirely real. (Note, for A 2 — | ^ the wavefunctions 
themselves can be complex, owing to the singularity of the potential at the origin and the 
departure of the quantisation contour from the real axis there.) As M moves away from 1, 
pairs of eigenvalues can become complex; as discussed earlier, this is always preceded by the 
coincidence of two real eigenvalues and so the first complex eigenvalues will emerge from points 
in the (2A, a) plane at which the spectrum has degeneracies for M = 1. We aim to investigate 
exactly how this occurs. 

In [27], Bender et al. used a perturbative approach to study the spectrum for M near 1 with 
a = and A 2 = \. The full Hilbert space was truncated to the subspace spanned by M = 1 
eigenfunctions |2n— 1) and |2n), where H\\m) = (2m+l)|m), m £ Z + , and Hm expanded within 
that two-dimensional subspace about Hm=i- Diagonalising the resulting 2x2 matrix yielded 
an approximation to the eigenvalues of Hm- However, as shown in [28], this approximation 
predicts level- merging for both signs of M—l rather than the one sign actually observed, and 
when applied to the pair of levels |2n) and |2n + 1), it predicts that they too will merge, contrary 
to the actual behaviour of the model. These problems can be traced to the fact that the M = 1 
eigenvalues at a = A 2 — \ = are equally spaced, making the truncation to the subspace spanned 
by |2n— 1) and |2n) unjustified. 

For the more general Hamiltonian (|5.2p the situation can be improved, as a and A can be 
tuned so as to make some pairs of levels close to each other relative to all of the others. Truncation 
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to these levels will then be reliable, and as we show below it gives a good approximation to their 
behaviour for M close to 1. 

To see how a consistent prediction of exceptional points can emerge from this approach, it 
is worth examining a simple 2x2 example which illustrates the main features. Consider the 
'unperturbed' Hamiltonian 

*(,) - (? 4) (5.3) 

where rj will be considered small but fixed, with the eigenvalues ±2?7 corresponding to the nearby 
pair of energies in the full problem. Add to it a perturbation with both diagonal and off-diagonal 
parts: 

e I a i \ 

VM = — . (5-4) 
■q \i -a I 

where a is fixed and e is the perturbing parameter (corresponding to M — 1 in the full problem). 
The factor of 1/rj will reflect the fact that nearby levels in the unperturbed problem interact 
more strongly as they approach each other. Then H\ +fi = H\ + V e has eigenvalues 



E± = ±y/(2r] - ae/rj) 2 - e 2 /r, 2 (5.5) 

and exceptional points at e = ±j^ry 2 . For fixed a / il the two exceptional points are at 
e = 0(r/ 2 ), so, even with the 1/rj factor in its specification, V e (if) is still small at their locations. 
For a = ±1 one exceptional point is pushed away to infinity, but the other remains in a region 
where the perturbation is still small. 

5.2 Perturbative locations of the exceptional points 

Returning to the original problem, the Hamiltonian at M = 1 is 

Hi = ~^ + x 2 + ^^-a. (5.6) 
With the given boundary conditions, H\ has c-normalised eigenfunctions [29] 

/n =S= = sV^ e -£ L ±V), ra = 0,l,... (5.7) 

a/(1 - e^ 2mX )T(±X + n + 1) 

where the Ln are Laguerre polynomials. The corresponding eigenvalues are 

£± = -a + 4n + 2 ± 2A . (5.8) 

A degenerate eigenvalue occurs when = for some n and m, which requires 

A = m — n . (5-9) 

Thus, on the vertical lines 2A € 2Z in the (2A, a) plane, infinitely-many pairs of the eigenfunc- 
tions (|5.7p are proportional to each other. Indeed, if A = q is a non-negative integer, then for all 
non-negative integers p, 4>p +q = i(—l) g (pp. Since cp^ — ► <fi~ when A — > —A, it also follows that 
4>p+ q oc 4>p when A = —q. 

In order to find the eigenvalues of Hm for M = 1 + e, we treat Hm = H\ +e in a basis of 
near-degenerate eigenfunctions of H\ by writing it as 

H 1+e = H 1 + V e (5.10) 

where Hi is given by (15. 6p and 

K = -x 2 - (ix) 2+2e - a{ixf. (5.11) 
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The exact matrix elements of V e in the truncated basis of H% eigenfunctions were found by 
Millican-Slater [29], and are reproduced in appendix [Cj while those of Hi are given by (|5.8p . 
Rediagonalising the resulting 2x2 matrix gives the approximate energy levels. 

To find the exceptional points reliably, we require both that the perturbation is small, and 
that the two levels in the truncated subspace are close. With M = 1 + e and A = q + ry, this 
means that e and rj must be small. In fact, we shall see that the exceptional points occur when 
e is of order rj 2 , and our approximations will be good in this region. We shall also assume that 
q > 0, as results for negative q are easily restored using the A — > —A symmetry of the problem. 
For small values of rj, the pairs of levels , (j)p+q}j P > 0, are almost degenerate; to lighten the 
notation, we fix the integer p > and denote the corresponding basis by {4> + , 4>~} = {4>p, <t>p+ q }- 
The matrix elements of H\ are 



\Hi\ 
\Hi\ 

m 



4p + 2q + 2r] + 2 
4p + 2q-2r] + 2 
<<n#l|0 + > = , 



(5.12) 
(5.13) 
(5.14) 



while those of V e follow from (j02]) . (IC3|) and (jOHJ) . 

Expanding in e and r\ and retaining terms proportional to rj, e/r], e and e 2 /r/ the matrix 



elements H. 



ah 



\HmW) are 



H 4 



4p + 2q + 2-a+(2p + q+l--)- + 27] 

\ 2/ rj 

+ ^2p + g+l-|)v(P + Q + l) + 2p + 2)e 

2 

+ ((2p + q + 1 - |) i>{p + q + 1) + 2p + l) - ; 



(5.15) 



4p + 2q + 2 
+ 



a 



2p + q+l 



n 



2r/ 



2p + g + 1 - - ) ij){p + 1) + 2p + 2g + 2 ) e 



(5.16) 



n 



2p + g+ l- - r )V(p + 9 + l) + 2p + l) — ; 



2 

H 
+ 

+ 



2p + 9 + 1 - - ) e 



-{2p + q + l-- 



) + 9 + 1) - ip(p + 1)) - Q J er/ 



(5.17) 



2p + q + l 



a 



ip(p + q+l) + 2p + l)e 



where ip(z) = T'(z)/T(z). Diagonalising H^, the approximate eigenvalues E± at M = 1 + e, 
\ = q + rj, a are: 

#± = 4p + 2g + 2-o+ (2p + q + 2 + ^(4p + 2q + 2-a)(ij(p + q +l)+ip(p+l))je 



± 



(8p + 4g + 4-2a)e + 4r/ 2 

+ ((4p + 2g + 2 - a)(i/)(p + q + l)-ip{p+ 1)) - 4g) er] 

1/2 

+ 4<? + 4 - a)^(p + q + 1) + 8p + 4 ) e 2 



(5.18) 



Within this approximation, exceptional points occur on the curves on the (2A, a) plane where the 
argument of the square root in (|5.18p vanishes. These curves, and their images under A — ► —A, 
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are plotted in figures [19] and [20] for e = 0.005 and e = 0.02 respectively. Each shows the excep- 
tional lines corresponding to p and q equal to 0, 1 and 2. (The exceptional lines for other values 
of p and q are outside the regions shown on the plots.) The dotted lines indicate a± S Z + , as 
previously. 
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Figure [T9l Perturbative lines of exceptional points for M = 1.005. 
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Figure [201 Perturbative lines of exceptional points for M = 1.02. 

As M increases, regions of complex eigenvalues open up from the lines A £ Z, starting near 
the bottom of the spectrum. While the mergings of these regions and the joinings of their 
exceptional lines to form cusps cannot be seen within this approximation (since the truncation 
is to just two levels), the pictures are consistent with the numerical evidence in the last section 
that the cusps move down from a = +oo as M increases from 1 towards 3. 
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The clearest insight into the transitions near M = 1 comes on retaining only the leading 
terms of the matrix elements for small r\ and e, namely those proportional to r/ and e/rj. For 
A = q + rj and M = 1 + e as before, the matrix elements in the basis {4> + , (/>"} = {<t>p, ^p+q} 
simplify to 






-ne/rj —ine/r] 
-ine/r] ne/rj 



where 



The approximate eigenvalues are then 



2 a 



2p — q — 1 . 



E, 



approx 



-2n±2y/r] 2 - ne 



(5.19) 



(5.20) 



(5.21) 



Apart from the overall shift by —2k and the replacement of e by ke, (15. 19ft and (15.21H have 
exactly the same form as the toy example (|5.5p at a = 1, one of the two values for which an 
exceptional point is found for only one sign of e. Thus our approximation captures an important 
feature of the full problem which was missed by the simpler approach used in [27] . Exceptional 
points occur when the argument of the square root in (|5.21j) vanishes. At fixed e, and using the 
A — > —A symmetry, this happens on the parabolas 



1 

a = 4p + 2q + 2 + — (2A ± 2q)' 



(5.22) 



on the (2A, a) plane, where p and q are non-negative integers. Thus there is a parabola rooted 
at every intersection of the lines a + £ Z + , a~ € Z + . However, there is a significant difference 
between the situations for e > (M > 1) and for e < (M < 1). For e > 0, the parabolas are 
upwards convex, as in figures [T9l and [20l above. Any fixed value of A and a in the neighbourhood 
of a line A = q within which the 2x2 truncations are valid lies inside only finitely many of the 
parabolas centred on that line, and thus sees only finitely many complex eigenvalues. 



a 




Figure [21] Perturbative lines of exceptional points for M — 0.98, 
with only a subset of the lines shown. 

By contrast, for e < the parabolas are oppositely-oriented, as in figure [2TJ Any given 
point (2A, a) near to a line A = q now lies inside infinitely many of the parabolas centred on 
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that line, and outside only a finite number of them. Thus truncation predicts that infinitely- 
many eigenvalues will be complex, with only finitely many remaining real, these real levels 
lying at the bottom of the spectrum. This is exactly as is observed in the full problem. The 
transition to infinitely-many complex eigenvalues was first noted by Bender and Boettcher [3] 
for A = 1/2, a = 0. It was subsequently treated analytically, for general A though still with 
a = 0, in [28], using a non-linear integral equation for the eigenvalues found via the Bethe 
Ansatz approach to the problem. While the latter approach is more systematic, the perturbative 
understanding of the phase transition just given is particularly transparent, and gives a more 
immediate understanding of the regions in the (2A, a) plane where complex levels are first to be 
found. 

A check on the truncation method can be made using the asymptotic obtained in [28] for 
the value of M = M cr i t < 1 at which high-lying eigenvalues E merge. With M cr { t = 1 + e cr ; t and 
A / 1/2, this iaj: 

41n|cos(7rA)| , . 

e C rit ~ • (5.23) 

For X = q + 7] and r\ small, this implies e cr it ~ —2i] 2 /E. This is easily seen to match the result 
just obtained, since (|5.21|) places the exceptional points at e = rf / 'k, and for E large, E ~ —2k. 

In table [2] the various approximations used in this section are compared with numerical data 
obtained from a direct solution of the ordinary differential equation. The numerical eigenvalues 
found by solving the full problem are denoted by E exac t; their numerical errors are smaller than 
the last quoted digit. The result using the 2x2 truncation and the exact matrix elements is 
Etrunci the initial approximated truncation (including the terms proportional to e and e 2 /r/) is 
E±, and the final approximation (retaining only terms proportional to 77 and e/rj in the matrix 
elements) is E approx . The table shows the comparison for sample values of e, a and 77, for for 
p = q = (i.e. A = 77) and p = 0, q = 1 (i.e. A = 1 + 77). 



p = q = 




e = 0.001, a = 0.9, rj = 0.01 


e = 0.001, a = 0.9, 77 = 0.25 


Eexact 
Etrunc 

E± 

Eapprox 


1.05069482 
1.05069431 
1.05067066 
1.04900980 


1.15266823 
1.15266441 
1.15269439 
1.15099019 


0.599733995 
0.599733083 
0.599485149 
0.597804819 


1.60332810 
1.60332370 
1.60387991 
1.60219518 


p = 0, q = 1 




e = 0.01, a = 3.9, 77 = 0.01 


e = 0.01, a = 3.9, 77 = 0.25 


Eexact 
Etrunc 

E± 

E 

^approx 


0.07899348 
0.07897778 
0.07913480 
0.05101020 


0.18089945 
0.18034086 
0.18078797 
0.14898979 


-0.36215520 
-0.36225580 
-0.36280969 
-0.40199601 


0.62170890 
0.62111404 
0.62273248 
0.60199601 



Table 2: Comparison of the various approximation methods used for MkI. 



The treatment so far has concerned the limiting region |e| <C 77 <C 1, which suffices to capture 
the behaviour of the exceptional lines as 77 — ► 0. Other limits are also interesting, and in closing 
this section we remark that other presentations of the Hamiltonian may then be useful. As an 
example, we return to the toy model (|5.3p . (|5.4p . at a = 1, and consider taking 77 — » before 

"•"When comparing with eq.(5.37) of [28], note that the e used there is equal to 2M — 2, and not M — 1. 
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e — * 0. As in [17], one can introduce a pair of matrices 



V2\i 1 



i? 



9 



(5.24) 



where q 2 = 2ie/rj. Then H\ +e = H\ + V t is similar to 



/ l + rf/e 
+ + Ue 



(5.25) 



It is now possible to set rj = 0, showing that the Jordon block is indeed recovered as the limit is 
taken. 

6 Perturbation theory about M = oo 

In this section we complete our analysis with a perturbative study about the model at M = oo, 
which is shown in appendixOto be exactly solvable. For large M, the second duality of appendix 
FPU maps the original eigenproblem (jl.ip into the Schrodinger equation 



where 



M 



-1 + 



M + 1 



d2 ±, \ 



1 + e , E 



z 2 + 



X 2 



a 



s 2M 

/ 2 \w 
\M+1 J 



E 



1 



z) = -^E(-iz) 2e <j)(z) 



A 



Af+1 



A 



(6.1) 



Af+1 



a 



(6.2) 



and we have set e = 2/(M + 1). Under the duality transformation, the contour C transforms 
into a curve equivalent to an M-independent straight line running just below the real axis. 

The inhomogeneous complex square well of appendix |D] appears from (16. lj) in the large-M 
(small e) limit, when the right-hand side reduces to an additional angular momentum term so 
that (|6.ip becomes the (PT-symmetric) simple harmonic oscillator, when viewed as an eigen- 
problem for q. The (unnormalised) eigenfunctions 



correspond to the a eigenvalues 



A=\^X 2 + E n 



(6.3) 



c£ = 4n + 2±2A/A 2 + £„ 



(6.4) 



Alternatively the problem at M = oo can be considered at fixed a as a generalised eigenproblem 
for E, with the (entirely real) spectrum following on rearranging (|6.4p : 



E n = (2n + 1 



\ 2 + n = 0,1,... . 



(6.5) 



The pair of levels E n and E m , n ^ m, will be degenerate whenever d = 2(n + m+l). Thus 
degeneracies occur in the spectrum on the horizontal lines a = 4, 6, 8, . . . in the (2A, a) plane, 
and a perturbative treatment will be reliable close to these lines. 

The eigenvalue problem at large but finite M can be explored by taking e small and truncating 
the full Hamiltonian H e to the 2x2 subspace spanned by the eigenfunctions ^ associated with 
the levels 



(q-2p + r]/2) 2 -A 2 



(q-2p- r?/2) 2 - A 2 



q G Z H 



(6.6) 
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This pair of eigenvalues will be almost-degenerate when a = 2(q+l)+rj andp = 0, 1, . . . [(q—l)/2] 
provided rj is small. When rj is zero the eigenvalues merge to the single eigenvalue Eq := = 
E~ p . Since the eigenfunctions (|6.3p satisfy the nonstandard eigenproblem 

Hot* = —E^^ , (6.7) 
z z 

the usual inner product must be weighted by a factor of z~ 2 , and so we define 

4> n (z)(f) m (z)z~ 2 dz (6.8) 



with a small positive s to avoid any singularities at z = 0. Using the integral (IC.lj) and analytic 
continuation as necessary, the orthonormal eigenfunctions are 

<f>+(z) = V / 2p!(g-2p + r ? /2) ^-2^/2-^/2^2^/2^ (6 Q) 

^/(l-e^jr^-p + r/^ + l) P 

and 



dr(z) - 6+(z)\ - V 2 (l-P)K2p-<l + ll/2) l/2+2p- q+V /2 e -zV2 L 2p- q +v/2 ( 2) f g w) 

In the truncated basis any eigenfunction 4> can be approximated as <j) = fi(p + + v$r for some 
constants [i and v. Applying H e to <j), the corresponding approximate eigenvalue E must satisfy 



+ v E-^ = —(-i z )^E(p(P + + (6.11) 



1 

z* z n z 1 

given that (j)^ are eigenfunctions of the unperturbed Hamiltonian (|6,7p . Thus taking inner 
product of (16. lip with (j) in turn, we obtain 



2c 



E ( )I-!)_„-j2eU+( ;j2eu-l ) ( „ ) ' (6.12) 



We use the integral (jC.lf) in appendixOto calculate the required exact matrix elements. To 
leading order in e and rj, the matrix elements H a b = (r/> a |(— iz) 2e \(j) b ) are 

H++ « 1 + — + (f/>(g-2p) + ^(c7-2p+l) - e 

+ 2(V(g-2p) + V(g-2p+i) - ^(?-p+i)) — ; (6.13) 

H — « 1- — + ^(«-2p)+^(g-2p+l)-'0(p+l))c 

- 2^(g-2p) + V(g-2p+l) - ^(?-p+l)) — ; (6.14) 



#4 



+ 2^(?-P+l) -ip(q-2p) - ^( g -2p+l)) — ] • (6.15) 
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Diagonalising the RHS of (|6.12p . the approximate eigenvalues at e = 2/(M + 1) and a 
2(q + 1) + 77 are 



E± 



L\, + ( ^ (Hq-p+l) + t/>(p+l) - 4^(g-2p+l)) - ^ - 2(g-2p) ) r 



± 



4(g-2p)£ e + (<?-2p) V + (<?-2p)(V>(<?-p+l) - V(p+l))-Eoe»7 



1/2 



+ ( 2(g-2p) (40(g-2p+l) - 3^(p+l) - + 40(g-2p)) £ + 4(g-2p) 2 J e 2 

(6.16) 

where Eq = (q — 2p) 2 — A 2 . Just as for M ~ 1, the exceptional points can be located by finding 
where the argument of the square root in (|6.16p vanishes. Figure [22] shows the resulting curves of 
exceptional points in the (2A, a) plane for M = 250, taking q = 1 . . . 5 and p = . . . [(q — l)/2] . 
The match with the results from a numerical solution to the full problem is excellent, and indeed 
even at M = 30 the truncation method gives a plot essentially indistinguishable from that shown 
earlier in figure [T8l 



a 



1500-f 



1000- 
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Figure [5U Perturbative predictions for the exceptional lines for M = 250. 

The main features of the transitions are most clearly understood if only the terms propor- 
tional to 77 and e/77 are kept in the matrix elements H a t>. Rediagonalising (|6.12p . the approximate 
eigenvalues are 

(6.17) 



E, 



approx 



£ ± Viz - 2p) V - 4E (q - 2p)e . 



Demanding once again that the argument of the square root vanishes leads to the prediction 
that the exceptional points lie on the ellipses 



a 



M + l 



1 ( g -2p)-4e [(q-2p) 2 -4 



"(M + l) 2 







(6.18) 



in the (2A, a) plane. Thus as M decreases from infinity isolated ellipses of unreality appear, 
starting from segments of the degenerate lines a = 4, 6, 8, . . . at M = 00 and acquiring exactly 
the 'nested' structure seen in figures [18] and | 
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Table [3] compares the various levels of approximation used in this section with numerical data 
obtained from a direct treatment of the ordinary differential equation, in the same notation as 
table [21 The table shows the comparison for sample values of e, A and rj, for p = 0, q = 1 (i.e. 
a = 4 + rj) and for p = 0, q = 2 (i.e. a = 6 + rj). 



p = 0, q = 1 




e = 0.001, A = -1.2, rj = 0.01 


e = 0.001, A = -1.2, rj = 0.25 


Inexact 
Etrunc 

E± 

E 


-0.48512051 
-0.48511885 
-0.48514998 
-0.48312771 


-0.3988914 
-0.3988926 
-0.3989179 
-0.3968723 


-0.679527538 
-0.679526948 
-0.695319170 
-0.693495562 


-0.17313365 
-0.17313370 
-0.18874878 
-0.18650444 


p=0,q = 2 




e = 0.01, A = —3, 77 = 0.01 


e = 0.01, A = -3, rj = 0.25 


Eexact 
Etrunc 

E± 

E 

'-'approx 


-5.59855357 
-5.59525382 
-5.60435172 
-5.63277168 


-4.36272053 
-4.36529048 
-4.35836985 
-4.36722832 


-5.72941123 
-5.72605257 
-5.75706543 
-5.80622578 


-4.20419136 
-4.20649571 
-4.20565613 
-4.19377423 



Table 3: Comparison of the various approximation methods used in section [5] 



7 Conclusions 

In this paper we have continued the project initiated in [1,2], and mapped out the phase diagram 
of a three-parameter family of PT-symmetric eigenvalue problems related to the Perk-Schultz 
models. Special features have enabled us to make precise the Jordan block structures at a 
subset of the exceptional points, going beyond the finite-dimensional examples which were the 
subject of most previous work. We have also uncovered some novel properties of the Bender- 
Dunne polynomials. The resulting phase diagrams at fixed M, consisting of lines of quadratic 
exceptional points punctuated by triply-exceptional (cubic) cusps, generalise the previously- 
observed story at M = 3 in an appealing way, and the perturbative treatment about M = 1 has 
allowed us to understand the transition to infinitely-many complex eigenvalues which occurs as 
M decreases below 1 from a new perspective. The dualities that we have used were crucial in 
making a reliable numerical treatment of the problem, and may be of independent theoretical 
interest, especially given the roles that this set of models plays as possibly the simplest example 
of an ODE/IM correspondence. 
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A Basis for an n x n Jordan block 



(See [16] for a discussion of the n = 2 case.) To illustrate a method we can use to construct the 
basis of an n x n Jordan block, which arises when n eigenstates merge, we will work with a toy 
model. Take an n x n matrix L, depending on one parameter e: 



L(e) 

This has n independent eigenvectors: 

/ 



/0 1 

... 
\e 








... \ 



1 

o / 



(A.l) 



fa 



e (2Trij/n) e l/n 
e (Anij/n) e 2/n 

L e 2n(n~l)ij /n e (n-l)/n j 



n . 



(A.2) 



When e = 0, L(e) has a Jordan block form, but at this point all n eigenvectors tpj become equal 
and so no longer form a basis. We therefore need to construct a new basis consisting of the 
vectors <j>^ , k = . . . n — 1 which satisfy a Jordan chain 



L(e)0(°) 

For simplicity we begin with the eigenvector tp n where 

L(e)Vae) = e V >n(e) 



e=0 
e=0 



n 



e=0 



(A.3) 
(A.4) 

(A.5) 



Clearly <fr(°) = ipn{c) satisfies the condition ()A.3P when e = 0. We could choose <j>(°) to be any of 

(A.6) 
(A.7) 



the ipj here; each one would lead to a different normalisation for the <p > below. 
Before we construct the other basis vectors, we introduce some notation. Let 

n^l d 

D — ne n — 

de 



and 



L = 



dL 

~dl' 



Note that L is linear in e so 4^ = 0. We now have the following commutation relations 



[D,e h ' n ] 
[D,L] 
[D,L] 



ke (k-l)/n 

ne (n " 1)/n L 
. 



(A.8) 
(A.9) 
(A.10) 



Finally, define 



k + 1 



= (A.ll) 

By induction, it is easy to show that acting with D on (|A.5[) times for 1 < < n — 1 gives 
fc-i j 



En 

jf=0 i=0 



U + 1) 



-i-i) + £0(*)=0(*-i) + ^0(*). 



(A.12) 
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When e = this satisfies (|A.4|) . so an appropriate basis is 
and 



for k = 1 . . . n — 1. 



(fe) = I^-i) 



(A.13) 
(A.14) 



B Two dualities 

As noted in [30], useful relations between spectral problems which arise in the ODE/IM corre- 
spondence can often be found by simple variable changes. Here, starting from (II. ip and setting 
z = ix as in (|2.14p to obtain 



dz 2 



1>(*) + 



z lM + az M-l + ^_l + E 







(B.l) 



we exploit the fact that, for arbitrary j3, the combined substitutions z = y@, ifi(z) = y^ 3 1 " 2 0(3/) ) 
transform d 2 ijj/dz 2 without introducing a first derivative term: 



3/2-3/3/2 



d k \ y 



dy 2 



P 2 



Ay 2 



so that the equation becomes 



y 2(M+l)/3-2 + (Af+l)/3-2 + P X 



2 \2 _ 1 

4 i zp„ ,2/3-2 



(3 2 y 2 



+ Ey z 



0(y) = o 



(B.2) 



(B.3) 



Two important special cases are (3 = 1/(M+1) and (3 = 2/(M+l). 

1) = 1/(M+1) : setting y = kw with n = ((M+1)/\^E) M+1 leads to 

d 2 



dw 2 



{w) + 



w 2 



+ E 



w) =0 



where 



M 



M 



E 



(M+l 



A , a 



M+l' (-E) M + l ' M+l'"' ~ M+l 

This generalises the duality used in [30] to inhomogeneous potential^]. 
2) (3 = 2/(M+l) : setting y = kw with k = y/(M+l)/2 yields 



a . 



(B.4) 



(B.5) 



where 



M 



-1 + 



dw 2 



M+l 



(w) + 



, E 



w 2 + Ew 2M + 



+ a 



w- 



(w) = 



2M 
2 \ M +! 



m+i; 



E , A 



M+l 



A , a 



M+l 



■ a . 



(B.6) 



(B.7) 



To obtain an equivalence between eigenvalue problems, the transformation of the boundary 
conditions under the mappings must be tracked. The boundary conditions from section [1] trans- 
late into the requirement that eigenfunctions of the initial problem (IB.lj) should decay in iS—i 



§It is interesting that, while [30] is indeed the first time that this duality was applied in the context of integrable 
quantum field theory, the homogeneous case can be traced back to (Isaac) Newton: see [31,32]. 
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and iSi, where the sectors S k were defined in (jl.3p . After the transformation the simultaneous 
decay should instead be in iS-\ and iS\, where for case 1, Newton's duality, 

S k = {x e C : |arg(ia?) - vr/c | < vr/2 } , (B.8) 

while for case 2(0 = 2/(M+l) ), 

S k = {x G C : \&rg(ix) -irk/2 \ < vr/4} . (B.9) 

In both cases the transformed sectors are independent of M, reflecting the fact that the leading 
terms in (lB.5h and (jB.6[) at large \w\, E and w 2 respectively, are themselves independent of M. 
For the first duality it might appear that the sectors iS±\ coincide, but this is not so - the 
branch cut in the original problem (11.11) becomes a cut along the negative real axis of the w 
plane, and so the two sectors lie on top of each other on the full Riemann surface of the problem. 
For the second duality the sectors are those of the simple harmonic oscillator and this makes 
(|B.6P particularly useful for numerical work: eigenvalues can be found by solving the ODE on 
a straight, M-independent contour, running vertically (parallel to the imaginary axis) in the 
right half of the complex w plane. An efficient approach uses WKB asymptotics at large \w\ as 
initial conditions for a pair of numerical solutions, <f)-i and <pi, decaying as Qmw — > ±oo, and 
then locates the eigenvalues by looking for zeros of the Wronskian W[</>_i, (pi], evaluated in the 
neighbourhood of the origin where both numerical solutions are reliable. This method was used 
to produce many of the figures in this paper. 

Replacing why w/i trivially rotates the dual problems back to a more usual "PT-symmetric' 
form. The mappings can also be used to give equivalences for spectral problems initially specified 
by the simultaneous decay of eigenfunctions on more widely-separated pairs of Stokes sectors 
than S-i and S\. The homogeneous cases of these problems were discussed in [27], and related 
to fused transfer matrices in integrable models in [11]. 



C Useful formulae 

This appendix records a number of formulae used in the main text. All can be inferred from 
the following basic integral, involving a pair of Laguerre polynomials: 

e t i.i+p)^ e -t L P m {t)Ll{t)dt 

= — j — j r(i(7+p)+l+a) x 

n\ ml z 

3 F 2 (-m, \(p+^)+l+a, ±(p-j)+l+a; p+l, l(p-j)+l+a-n; 1) 

r(i( 7 +p)+l+q) ^ 
ml n\ 
m f \ 

( ? T)(p+l+A:) m -fe(|(p+7)+l+a)fc(^(/5-7)+l+a)fe(i(7-/o)-a)n-fc (G.l) 

where (a) n = a(a+l) . . . (a+n— 1) is the Pochhammer symbol and 3F2 is a generalised hyperge- 
ometric function. The first version of this result can be found in [29]; it generalises a formula 
for the case 7 = p that was given in [33]. The symmetry of the final expressions under the 
simultaneous exchanges m <-» n, p *-* 7 is not obvious, though it can be checked. 

In section [5] the matrix elements (4>^(x)\ (ix) 2M |<^(x)) and (cj)^(x)\(ix) 2M \(j)^ l (x)) were 
needed for general M, where <fin( x ) an d 4>n( x ) are the normalised wavefunctions given by (|5.7p . 
The relevant calculations were also carried out by Millican- Slater in [29], and we reproduce his 
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final results here. The matrix element (cf)^\(ix) is 
(cos(Mvr) + sin(M7r) cot(Avr)) 



-Af)„(A + l) 



v / n!m!r(A + m + l)T(X + n + 1) 
r(A+M+l) 3 F 2 (-m, A+M+l, 1+M ; A+l, 1+M-n ; 1) . 



(C.2) 



For M = 1 (one of the cases needed) there is a negative integer in one of the second group 
of entries of the hypergeometric function in (|C.2|) . and so for certain values of n and m these 
functions may be undefined. This is the case when n — 2 < m. However, for \n — m\ > 2 the 
symmetry of the inner products in n and m can be used to avoid the problem. In these cases, 
when M = 1, the (-M) n in the expressions above become (— l) n — so the inner products are 
zero. For n = m and n = m ± 1, by taking the limit M — ► 1 in (|C.2|) it can be shown [29] that 
the only non-zero inner products are 



{<j>+{x)\x 2 \<j>+{x)) 
(cP+ +1 (x)\x 2 \cPt(x)) 



1 + A + 2n 



n + 1 
n + A + 1 



(A-n) 



(C.3) 
(C.4) 



The matrix elements corresponding to (IC.2j) . (|C.3h and (IC.4D for n can be found by sending 
A^-A. 

The matrix element (^(x)|(ix) 2M |(/>~ (x)) is given by 



x )\{ixf M \r m {x)) = l - 



sin(M7r)(l - A) m (A - M) n T(M + 1) 



| sin(7r A) | i/r(l - A + m)r(l + A + n)m!n! 
x 3F2 (-m, 1 + M,M + 1- A; l-A,M + l- A- n; 1) , 



(C.5) 



which, unlike (|C.2|) . is always well defined at M = 1. 



D The inhomogeneous complex square well 

In the main text, the large-M limit of the spectrum of 



dx 2 



(ixf M -a(ix) M - 1 + ^- 1 ^ 

x A 



t/j(x)=Eil>(x), 4>(x)eL 2 (C), 



(D.l) 



was needed. The a = 0, A = 1/2 case was investigated in [34], where it was dubbed the 'complex 
square well'. To treat the more general case, we start with the same variable change as in [34], 
and set 

x=(-i + ^)E2U. (D.2) 



Taking the limit M — ► oo, using the identity limjvf— >oo(l + x/M) M = e x and dropping all 
subleading terms, (ID.lh becomes 



^ + ^E(1 + e™) + -VEa e^ 2 + — A 2 
dz 2 16 V ; 16 16 



V>0) = 



where 



\M+l) 



and the scaled parameters 



2A 
M + 1 



o 



2a 
M + 1 



(D.3) 

(D.4) 
(D.5) 
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were used to ensure the survival of the inhomogeneous and angular-momentum terms in the 
limit. Notice that in terms of A and a, the parameters a± of (|2.1f) are simply 



a± 



1±2A 



(D.6) 



The special feature of this limit is that the resulting ODE (|D,3j) is exactly solvable. Here we 
highlight the link with the simple harmonic oscillator by making a further variable change to 
w = £' 1 / 4 e l7r2 / 4 and trading ip(w) for 4>{w) = i/wip(w) . Substituting in, (f>(w) satisfies 



d 2 



dw 2 



+ 



2 , E + \ 2 -\ 
ur H k - 



cf> = —a 4> ■ 



(D.7) 



Boundary conditions should be imposed on the asymptotic Stokes lines z = ±2 — iy, y — > 00 [34], 
which translate into the positive and negative imaginary axes in the complex w plane. That 
said, the spectrum of (|D.7h can be recognised as that of the 'PT-symmetric simple harmonic 
oscillator [11,26], with 'energy' a and 'angular momentum' — 1/2 ± yE + A 2 (the reversed sign 
of the energy is a result of the rotated quantisation contour for (|D.7|) compared to that used 
in [11] ). Hence, from [11], (|D.7|) has a wavefunction normalisable on the quantisation contour 
if and only if 

(D.8) 
00 limit: 



a 



4n + 2±2y£ + A 2 , n = 0,l,... 
which translates into our main result for the exact spectrum of (jD.ip in the M 



E n = (2n + 1 



) - A 



n 



0,1, 



(D.9) 



Via E n = 4-E n /(M+l) 2 , this result also gives the leading behaviour of the original levels E n as 
the linit is taken. For A = a = 0, this reproduces the result of [34]. Notice that the spectrum is 
entirely real for all values of A and d, matching the situation at M = 1, the other exactly-solvable 
point. 
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